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ffi ■ ABSTRACT 

^ f It is generally thought that the light coming from the inner plunging region of 

black hole accretion discs contributes negligibly to the disc's overall spectrum, i.e. the 
plunging fluid is swallowed by the black hole before it has time to radiate. In the 
standard disc model used to fit X-ray observations of accretion discs, the plunging 
region is assumed to be perfectly dark. However, numerical simulations that include 
the full physics of the magnetized flow predict that a small fraction of the disc's 
total luminosity emanates from the plunging region. We investigate the observational 
consequences of this neglected inner light. We compute radiative transfer based disc 
spectra that correspond to 3D general relativistic magnetohydrodynamic simulated 
discs (which produce light inside their plunging regions). In the context of black hole 
spin estimation, we find that the neglected inner light only has a modest effect (this 
bias is less than typical observational systematic errors). For rapidly spinning black 
holes, we find that the combined emission from the plunging region produces a weak 
power-law tail at high energies. This indicates that infalling matter is the origin for 
some of the 'coronal' emission observed in the thermal dominant and steep power-law 
states of X-ray binaries. 



Key words: accretion, accretion discs ~ black hole physics - MHD - methods: 
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INTRODUCTION and spinning BHs (MacFadyen & Woolsley 1999), models of 

quasi-periodic oscillations (Abramowicz & Kluzniak 2001), 



The ubiquity and simplicity of black holes in our universe 
make them truly marvelous objects of study. The com- 
plete physics of each and every black hole (BH) can be r , , , , . , , 

, , . , 1 ,,,,,, , relatively straiehttorward; so lone as one can obtam the 

distilled down to lust two numbers : the black hole s mass , , , /• . ,. ,. 



must connect in some way to these two numbers. 

For binary systems, the task of measuring mass is 



period, orbital velocity, orbital geometry (i.e. inclination 

and eccentricity), and mass for a companion orbiting a 

-2 ; \ rr^n ■ • 1- 1 i, , ■ i , ■ , black holc, one can immediately obtain the black hole 

J/{GM c). Inis implies that all physical theories that . , , ■ ■, ,r , , 

. , , , , , , , ,. , , T-,TT ■ n mass using only Newtonian gravity (tor a recent example, 



M and angular momentum J, the latter of which is usu- 
ally expressed as a dimensionless spin parameter a* = 



involve black holes, e.g. the link between BH spin and 
jet power (Blandford & Znajek 1977), gamma ray bursts 



see Orosz et al. 2011b). The game of measuring black hole 
masses has been played as early as 1972 for Cygnus X-1 
(Webster & Murdin 1972; Bolton 1972), and to date we have 
robust mass estimates for about 20 stellar mass binary black 
hole systems (Remillard & McClintock 2006; Orosz et al. 
2007, 2009, 2011a,b; Cantrell et al. 2010), and 64 supermas- 
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1 In principle, BH charge is an independent parameter; however, Despite this success in measuring black hole mass, spin 

we do not expect astrophysical black holes to retain any signifi- has been a more difficult quantity to obtain. The spin of 

cant charge. an object only makes itself known at very short distances 
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through a general relativistic effect known as frame drag- 
ging. Thus to probe black hole spin, we must rely on observa- 
tions of matter that is very close to the BH horizon. In prac- 
tice, the only available probe is the accretion disc orbiting 
the black hole, which transitions from nearly circular orbits 
to plunging at a special location known as the innermost sta- 
ble circular orbit (ISCO). By observing the light given off by 
the accretion disc, it is possible to determine this transition 
radius. When the ISCO is expressed in terms of the gravita- 
tional radius Vg — GM/(? , it has a well-defined monotonic 
dependence on black hole spin (e.g. Shapiro & Teukolsky 
1983). Therefore, a measurement of the ISCO size yields the 
BH spin if the mass is independently known. 

For a black hole binary system, we measure the size of 
the ISCO by modelling the light given off by the accretion 
disc. Objects inside the ISCO (hereafter referred to as the 
'plunging zone') cannot remain in stable circular orbits, and 
are forced to plunge into the black hole in just a few free-fall 
times. Since the plunging time-scale is short, it is thought 
that fluid inside the plunging zone does not have time to ra- 
diate (Page & Thorne 1974), which means that if we could 
resolve an image of the accretion disc, then we would see a 
dark void corresponding to the plunging region in the cen- 
tre. In practice, we cannot resolve images^ of accretion discs 
around black holes, and thus our only information comes in 
the form of X-ray spectra. 

One method of estimating BH spin"* works by fitting 
the spectral shape of the thermal continuum emission from 
the accretion disc to thereby estimate the radius of the ISCO 
(this method was pioneered by Zhang et al. 1997a; see Table 
3 for some recent results). 

This 'continuum fitting' method uses the colour tem- 
perature of the thermal component, the distance to the 
source, and the received X-ray flux to obtain a character- 
istic emitting area for the disc. Given a model for the radial 
dependence of the disc emission, this characteristic emit- 
ting area determines the ISCO radius. The main drawback 
of the continuum fitting method is that we first need ac- 
curate estimates of the BH distance (to turn fluxes into 
areas), disc inclination (to turn the projected area into an 
ISCO radius), and BH mass (to get spin from the monotonic 
risco/vg vs. a* relation). Luckily, the techniques needed for 
measuring distance (Reid et al. 2011), mass, and inclination 
(Orosz et al. 2011b; Cantrell et al. 2010) for binary systems 
are well developed, and to date have been successfully ap- 
plied to more than half a dozen black hole binary systems 
(McClintock et al. 2011). 

The continuum fitting method has evolved through a 
sequence of progressively more complex disc models. The 
simplest model assumes multitemperature blackbody emis- 
sion from the disc (Mitsuda et al. 1984). Including relativis- 
tic effects such as Doppler beaming, gravitational redshift- 
ing, and light bending yielded the next generation of mod- 
els (called KERRBB - Li et al. 2005). The current gener- 

Although it has been recently proposed to use radio VLBI ob- 
servations to test GR by resolving the apparent shape of SgrA* 
(Doeleman et al. 2009). 

^ Another commonly applied method is known as the 'iron line' 
method, which estimates BH spin through spectral modelling of 
the Fc Ka fluorescence line (Fabian et al. 1989, see Miller 2007 
for a recent review). 
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Figure 1. Comparison of dimensionless luminosity profiles from 
the GRMHD simulations of Penna et. al 2011 (solid) with 
the standard disc model (dotted - Novikov & Thorne 1973; 
Page & Thorne 1974). Results from models A,B,C (defined in Ta- 
ble 1), are shown. The dashed line depicts an alternate measure 
of the GRMHD luminosity (described in §7.2 and Appendix D). 
Note that the luminosity in the standard NT disc model goes to 
zero at the ISCO. 

ation of disc models (named bhspec - Davis et al. 2005; 
Davis & Hubeny 2006a) frees itself from the blackbody as- 
sumption, and instead obtains the disc spectra by means of 
radiative transfer calculations. All of the above mentioned 
models use the prescription of Novikov & Thorne (1973, 
hereafter NT) as the underlying disc model'*, which assumes 
a perfectly dark void inside the ISCO. The advance that 
we make in this work is to apply radiative transfer calcu- 
lations to compute the spectrum for a general relativistic 
magneto-hydrodynamic (GRMHD) simulated disc that pro- 
duces emission from inside the ISCO. 

Currently, a crucial assumption in the continuum fit- 
ting enterprise is the perfect darkness of the plunging zone. 
This scenario applies only to the idealized case of a razor 
thin unmagnetized disc (see NT). For finite thickness discs, 
it has been argued that as long as the disc remains geomet- 
rically thin (with disc aspect ratios h/r <^ 1), the disc can 
be well approximated by the NT model (Paczynski 2000; 
Afshordi & Paczyrisky 2003; Shafee et al. 2008a). However, 
recent work with magnetized discs (Shafee et al. 2008b; 
Noble et al. 2009, 2010, 2011; Penna et al. 2010) suggests 
that there may be departures from NT even in the limit of 
thin discs. An observational difference between these mag- 
netized discs and the classic NT discs is the appearance of 
non-negligible luminosity inside the plunging region (see Fig. 

* This is the relativistic generalization of the standard 
Shakura & Sunyaev (1973) model for viscous geometrically thin 
but optically thick accretion discs. 
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1 for the result from Penna et al. (2010), where the extra 
plunging region luminosity constitutes ~ 4 per cent of the 
total). Recent models of unmagnetized discs that include the 
physics of energy advection (Abramowicz et al. 1988; most 
recently S§,dowski 2009, S^dowski et al. 2011) also share the 
feature of nonzero luminosity inside the ISCO, but in these 
advective models, the extra light only becomes significant 
when the accretion rate approaches the Eddington rate. 

The primary goal of the present project is to investigate 
the importance of the neglected light from the inner plunging 
region of accretion discs. Since GRMHD discs exhibit the 
phenomenon that we wish to investigate (i.e. nonzero light 
from the plunging region), we use them as the basis for all 
our investigations. We analyse the GRMHD simulations of 
Penna et al. (2010) in this work. To connect the GRMHD 
simulations with observables, we compute realistic (i.e. ra- 
diative transfer based) disc spectra from them. Essentially, 
we adapt the ideas pioneered in bhspec (Davis et al. 2005; 
Davis & Hubeny 2006a) to the domain of GRMHD discs. 
One limitation of the GRMHD simulations is that they do 
not inherently model the radiation field of the disc fluid. The 
way we put the photons back in (for the purposes of comput- 
ing disc spectra) is to do a radiative transfer post-processing 
step. To get the integrated disc spectrum corresponding to a 
particular simulation run, we take the following three steps: 

1) We slice the GRMHD disc into many individual annuli; 

2) For each annulus, we compute the local spectrum through 
a radiative transfer calculation; 3) By means of ray tracing, 
we sum up the light from every annulus to get the overall 
spectrum of the accretion disc. 

Several previous studies have also used GRMHD simu- 
lations to consider the impact of emission from the plung- 
ing region on disc spectra (e.g. Beckwith et al. 2008a; 
Noble et al. 2009, 2011; Kulkarni et al. 2011), but all as- 
sumed (modified) blackbody emission. Since the plunging 
region is precisely where the fluid becomes tenuous and 
optically thin, here the blackbody assumption is likely to 
be strongly violated. Our work extends the analysis of 
Kulkarni et al. (2011, hereafter K2011), which also made use 
of the GRMHD discs from Penna et al. (2010). The advance 
that we make beyond K2011 is that we obtain disc spectra 
by means of radiative transfer calculations (this allows our 
work to be free from the blackbody assumption which under- 
pins all previous work on GRMHD disc spectral modelling). 

Our paper is organized as follows: In §2 we describe 
the setup and properties of the simulated GRMHD discs 
that we use to generate disc spectra. We then detail the 
physics of the radiative transfer calculation in §3, and list 
the simulation quantities needed to generate the spectrum 
for each annulus. In §4, we explain how these simulation 
quantities are extracted, and we give a brief overview of 
the ray tracing process in §5. The implications from the 
additional plunging region light are discussed in §6, and in §7 
we briefly touch upon the limitations of our disc model. We 
summarize the key results in §8. Appendices A-D expand 
upon the technical details of various aspects of this work. 



HARM (Gammie et al. 2003; McKinney 2006; 
McKinney & Blandford 2009). The code works in di- 
mensionless units (G — c = ks = 1), and we assume 
that the fluid follows an ideal gas equation of state 
gas — (r-l)(7 where Pgas is the gas pressure, and U is the 
gas internal energy density. We choose the adiabatic index 
to be r = 4/3, which corresponds to an ultrarelativistic 
equation of state. We start off the simulation with a torus 
of gas that is initially in pure hydrodynamic equilibrium 
(De Villiers et al. 2003; Gammie et al. 2003) threaded by a 
weak (100 < Pgas/-Pmag < 1000) 4-loop poloidal magnetic 
field structure (see Penna et al. 2010 for details on the field 
topology). The torus is also set up such that its orbital spin 
axis is aligned with the spin axis of the black hole. 

There are three primary degrees of freedom for the 
torus: 1) the initial gas entropy, 2) the initial magnetic field 
strength and geometry, and 3) the initial angular momentum 
profile. The choice of initial gas entropy controls the overall 
thickness of the disc, wheras the choice of initial magnetic 
field strength/geometry controls the strength of the turbu- 
lence (Beckwith et al. 2008b). The angular momentum pro- 
file sets the radial extent of the starting torus. We evolve the 
system in time via a conservative Godunov scheme, with an 
additional caveat that we cool the gas via the following pre- 
scription: 

dU _ ^^ \n{K/K,) ^ 
dr T cool 

where U is the internal energy of the gas, K = P/ 
is the gas entropy constant, Ki — Pi/f^ is the initial 
entropy constant, and we set TcooI = 27r/Slfc for the cooling 
time- scale (f^fc = yfcWp). After a gas element heats 
up from dissipation, energy is removed according to Eq. 1 
such that the gas returns to its initial entropy^ (this acts to 
preserve the initial aspect ratio of the disc). In the absence 
of a full radiative transfer calculation, the cooling function 
is a substitute for the radiative energy loss that we expect 
from a geometrically thin, optically thick disc. 

The initial conditions of the specific simulations that 
we consider are listed in Table 1. All models are run using 
a fixed spherical polar mesh with Nr = 256 (radial cells). 
No = 64 (poloidal cells), and A''^ = 32 (toroidal cells), except 
for model E, which has twice the number of poloidal and 
toroidal cells. The cell spacing is chosen such that resolution 
is concentrated near the horizon and the disc midplane (see 
Penna et al. 2010 for details). We find that the increase in 
resolution from model D to E has the effect of increasing 
a. Models A-C correspond to thin discs and most of our 
reported results focus on these models. Model F differs from 
the other runs in its initial magnetic field topology in that 
it uses a single poloidal loop configuration rather than the 
4-loop model adopted in the other runs. 



2 GRMHD SIMULATIONS 

The simulation numerically evolves the 3D GRMHD s The cooling function in Eq. 1 is not invoked when K < Ki (i.e. 

equations in the Kerr spacetime via the code the cooling prescription does not add any heat to the fluid) 
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Table 1. GRMHD disc parameters 



Model a* Target h/r L/L^^^ a Comment 



A 


0.9 


0.05 


0.35 


0.22 




B 


0.7 


0.05 


0.32 


0.10 




C 


0.0 


0.05 


0.37 


0.04 




D 


0.0 


0.1 


0.70 


0.04 




E 


0.0 


0.1 


0.71 


0.08 


High Res. 


F 


0.0 


0.05 


0.36 


0.03 


1-Loop 



Note - The BH mass was set to be M = IOM0 . L/I/Edd and a 
are derived quantities that are measured from the simulations 
(see §4). K2011 also used Models A, B, C, and F in their 
work. 



3 ANNULI SPECTRA 

We obtain the overall GRMHD accretion disc spectrum by 
summing up the hght contributions from all annuli that com- 
prise the disc. Most previous work on disc spectra assume 
a modified blackbody prescription for the local spectrum 
emitted by each annulus. The blackbody is modified such 
that the specific intensity is given by: 

/.(T) = r^B^T), (2) 

where Bv{T) is the standard Planck function and / is spec- 
tral hardening constant (often assumed to be 1.7 for X- 
ray binaries). Following Davis & Hubeny (2006a), we im- 
prove on this assumption by solving the radiative trans- 
fer equation. The emergent spectrum for each annulus 
is obtained through the stellar atmospheres code tlusty 
(Hubeny & Lanz 1995). The code simultaneously solves for 
the vertical structure of a plane-parallel atmosphere and 
its angle dependent radiation field, and is independent of 
any blackbody assumption. The atmosphere is modelled as 
a series of ID vertical cells that are in radiative, hydro- 
static, and statistical equilibrium (which allows for depar- 
tures from local thermodynamic equilibrium in the atomic 
populations). In each cell, we obtain: the spectrum as a func- 
tion of viewing angle Iv{9), density p, gas temperature T, 
vertical height above the midplane z, and the particle com- 
position (rie, rijjj , Tijjjj , etc.). 

We are primarily interested in each annulus' emergent 
spectrum (i.e. Iv(9) for the surface cell). To obtain a unique 
solution in the stellar atmospheres problem, we need to spec- 
ify the following 3 boundary conditions: 

(1) The radiative cooling fiux F = asB.T^g (where Tcfi is 
the annulus effective temperature). 

(2) The vertical tidal gravity experienced by the fiuid 
(parametrized by Q = O-l/z where g±^ is the vertical 
acceleration and 2 is the height above the midplane). 

(3) The column density to the midplane m = E/2 (where 
E is the total column density of the annulus). 

The problem of computing the disc spectrum simply 
becomes one of obtaining the radial profiles of Tesir), Q{r), 
and m(r) for the GRMHD disc. We then compute the spec- 
trum for each disc annulus with tlusty using the corre- 
sponding values of TcS, Q, and m. 



3.1 Assumptions in the TLUSTY model 

In our annulus problem, we make the assumption of uni- 
form viscous heating per unit mass, and we ignore the ef- 
fects of magnetic pressure support, fiuid convection, and 
heat conduction. We allow for deviations from LTE by ex- 
plicitly computing the ion populations for H, He, C, N, 
O, Ne, Mg, Si, S, Ar, Ca, and Fe assuming solar abun- 
dances. In the non-LTE calculations, only the lowest en- 
ergy level is considered for each ionization state, except for 
the cases of H and He"*", which are modelled with 9 and 
4 levels respectively. The treatment of bound-free transi- 
tions includes all outer shell ionization processes, coUisional 
excitations, and an approximate treatment for Auger (in- 
ner shell) processes. Free-bound processes are modelled in- 
cluding radiative, dielectric, and three-body recombination. 
Free-free transitions are modelled for all listed elements, 
whereas bound-bound transitions are not modelled. Comp- 
tonization is included through an angle-averaged Kompa- 
neets treatment (Hubeny et al. 2001). Finally, to reduce the 
number of independent vertical cells in the problem, we as- 
sume that the annulus is symmetric about the midplane. 

One drawback of the algorithm used in tlusty 
(lambda-iteration and complete linearization, see 
Hubeny & Lanz 1995 for details) is that tlusty re- 
quires an initial guess for the conditions in each vertical 
cell. The code often fails to converge when the initial 
guess is poor. However, if tlusty finds a converged 
solution, it is robust to the choice of initial guess. Often, 
manual intervention (in the form of providing better initial 
guesses) is needed to ensure convergence, which means 
that the process of generating annuli spectra cannot be 
completely automated. However, since only three parame- 
ters are necessary to uniquely specify an annulus, we have 
simply manually precomputed a grid of annuli spanning 
the full range of parameter space for the case of stellar 
mass black hole systems. The grid spans logjg Tcft G 
{5.5, 5.6, 7.3}, login ^ {-4.0, -3.9, . . . , 9.0}, and 
logio m G {0.5, 0.6, . . . , 2.8, 2.9, 3.0, 4.0, 5.0, 6.0}. We then 
interpolate on this grid to obtain any needed annuli spectra 
(see Appendix C for details on the interpolation process). 



4 SLICING THE GRMHD DISC INTO ANNULI 

Since we wish to describe the accretion disc from the 3D 
GRMHD simulations as a series of annuli, we take out the 
poloidal and toroidal structure through an averaging pro- 
cess. All quantities are first subject to azimuthal averag- 
ing (since we model the annuli as azimuthally symmetric 
structures), followed by vertical averaging (i.e. averaging 
over 9) weighted by rest mass density. For each annulus in 
the GRMHD simulation, we must identify Tctf{r), Q{r), and 
E(r) before we can call on tlusty to provide the local spec- 
trum. 

We obtain Tcff (r) directly from the simulation luminos- 
ity profile (§4.1), Q{r) from the Kerr metric and the fiuid 
velocities (§4.2), and E(r) from the fiuid velocities and ac- 
cretion rate (§4.3). Since the simulations are dimensionless, 
we must first choose a black hole mass M and a disc lumi- 
nosity I//I/Edd to dimensionalize the radial profiles of Teff(r), 
Q(r), and E(r). For all the simulation runs, we have chosen 
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a fiducial value of M = 10 M©, and we set L/L^dd so the 
disc thickness in the radiative model matches that of the 
GRMHD simulation (see §4.1.1). 

One complicating factor is that the simulated disc is 
only reliable for a short range of radii. Beyond a certain 
critical radius r > ru (where Tie stands for the radius of 
inflow equilibrium - see Penna et al. 2010 for a detailed dis- 
cussion), the simulation has not been run long enough for 
the fluid to reach its equilibrium conflguration. The rele- 
vant physical time-scale for reaching inflow equilibrium is 
the accretion time face = r/u^ , where is the fluid's radial 
velocity. The problem with fluid outside of Tie is that this 
fluid still has memory of its initial conditions (i.e. fluid out- 
side of Tie has tacc > tsim , where tsim is the total runtime 
of the simulation). Beyond the inflow equilibrium point, we 
instead turn to an analytic disc model that is matched to 
the simulation disc (we use a generalized NT model that 
allows for nonzero stresses at the ISCO, see Appendices A 
and B for details). The analytic model also requires us to 
specify a disc viscosity parameter a (defined as the ratio 
between the vertically integrated stress and pressure, see 
Novikov & Thorne 1973). Table 1 lists the values of L/Z/Edd 
and a corresponding to the different simulation runs, where 
both parameters are determined by matching the GRMHD 
disc at Tie to the analytic generalized NT disc (see §4.1.1 
and §4.3.1). 



4.1 Flux profile - T^sir) 

From the simulations, we have extracted in the Boyer- 
Lindquist {t, r, 9, (f>) frame, the fluid four- velocity 
{u*^ jU"^ ,u^), and luminosity Lbl, where the lumi- 
nosity is computed from the cooling function dU/dr (see 
Eq. 1) in the following fashion: 



LBL(r) 



dU 



/ —gdtdr' d6d4>, 



</) = 9=0 r'^TH t = *i 



(3) 

where ti and tf are the initial and flnal times considered 
for the time integral, rn is the black hole horizon, ut is the 
fluid's conserved speciflc energy, and \/—g is the determinant 
of the Kerr metric. We do the integration for only the bound 
fluid, and then obtain the comoving flux by transforming 
from the Boyer-Lindquist flux (see K2011): 



0-SBrcfl('") = Fcomir) 



FBL(r) 



1 



diBL(r) 



-4:Tvrut dr 



(4) 



However, since the simulations are converged only for a short 
radial extent, we extrapolate the flux proflle to large radii by 
matching the GRMHD flux with an analytic disc model at a 
matching radius Tie. The extrapolation model we use is the 
same as K2011, which is a generalized Page & Thorne (1974) 
model (see Appendix A) that does not assume zero torque 
at the disc inner boundary, and hence allows for nonzero 
ISCO luminosities. 

To determine the point at which we can no longer trust 
the simulations (in other words, to locate r^e), we look at the 
radial profile of the GRMHD accretion rate (Figure 2). Be- 
yond a certain radius, we see that the accretion rate begins 
to significantly deviate from a constant value (the devia- 
tion from a constant accretion rate only occurs in fluid that 
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Figure 2. The GRMHD simulation accretion rates normalized 
to M at the horizon are shown. We identify the region left of the 
vertical dashed line as converged, so the line marks the boundary 
Tie where we match the GRMHD disc to an analytic model. 



has not reached infiow equilibrium). We pick the matching 
point to be the outermost point where we have well behaved 
accretion, though the final luminosity profile is fairly insensi- 
tive^ to our choice of matching radius (i.e. after varying the 



matching location by ±20 per cent, 
changes by less than 10 per cent). 



the luminosity profile 



4.1.1 Disc thickness matching (L/LEdd) 

The GRMHD luminosity that we measure (as shown in Fig- 
ure 1) is actually in dimensionless units of L/Msim where 
Msim is the local accretion rate of the simulations (see Fig- 
ure 2). To calculate Tcs in dimensioned units of [K], we need 
to first determine the dimensioned luminosity corresponding 
to the GRMHD simulations. Unfortunately, the simulations 
do not include radiation physics, so the simulation fiuid does 
not have a set density scale, which means there is no direct 
way of measuring M and hence L. We thus resort to an indi- 
rect method for determining this accretion rate. We identify 
L/L-Edd of the simulation with the value that matches the 
disc thickness in both the GRMHD simulated disc and the 
TLUSTY annuli'^ (see Figure 3). However, since tlusty ig- 
nores magnetic support and the simulations ignore radiation 
pressure, the tlusty based and GRMHD thickness profiles 
have completely different shapes. We opt to match the thick- 
ness at the radius corresponding to the luminosity peak (see 
Figure 3). 

Although the luminosity profile drops off fairly rapidly 
inside the ISCO (see Fig. 1), we find that the disc's effective 
temperature remains roughly constant (Fig. 4). Therefore 
the falloff in disc luminosity when approaching the BH hori- 
zon is purely a geometric effect. Annuli near the horizon 



^ The insensitivity of the luminosity extrapolation to the match- 
ing radius r^e was demonstrated in appendix A of K2011. 
^ K2011 also applied this method to dimcnsionalizc their lumi- 
nosity profiles 
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Model C - (a* = 0) 




Figure 3. Shown are the h/r disc thickness profiles for model 
C, as computed by TLUSTY and our GRMHD simulation, where 
h = J p\z\dz/ J pdz. The black dotted lines depict the TLUSTY 
disc thicknesses for different choices of luminosity (the four sets of 
dotted lines correspond to L/I/Edd = 0.1,0.2,0.4,0.8). The black 
solid lino denotes the matched thickness profile with L/LEdd = 
0.37. The black dot represents the radius where the luminosity 
is greatest, which we have chosen to be the radius at which we 
match the TLUSTY and GRMHD thickness profiles. 




Figure 4. Here we plot the T^g flux profile, computed from the 
luminosity profiles of Figure 1, for the accretion rates listed in 
Table 1. 



have less surface area, and hence they contribute less to the 
total luminosity. 



4.2 Vertical gravity profile - Q{r) 

The vertical gravity parameter g±, which represents the 
tidal vertical acceleration in the disc, is obtained by eval- 
uating the 7?|.j: component of the Riemann curvature tensor 
in the comoving frame. A commonly used prescription for 
the vertical gravity is that of Riffert & Herold (1995), which 
assumes that the disc fluid moves along circular geodesies; 



this approach becomes invalid inside the plunging zone. We 
thus turn to the prescription of Abramowicz et al. (1997), 
which relaxes the circular orbit assumption: 

g^{r,z) = nlR,{r)- z, (5) 

where ~ {GAd/r'^)^^^ is the Keplerian orbital frequency, 
and the dimensionless relativistic factor Rz(r) is given by: 



+ a'i{ut - 1) 



(6) 



In Eq. 6 we use the four-velocity corresponding to the mid- 
plane fluid from the GRMHD simulations for r<rie (as 
determined from Figure 2). For r>rie, we use circular 
geodesies as the four- velocity in Eq. 6. We find that the 
two prescriptions for vertical gravity (Riffert & Iferold 1995; 
Abramowicz et al. 1997) give nearly identical results outside 
the ISCO. 

For convenience, we define the function Q{r) such that: 



^ Q{r) ■ 



(7) 



and we interpret Q{r) as the radial dependence of the ver- 
tical gravity. 

4.3 Column density profile - E(r ) 

Since the GRMHD simulation does not include radiation 
physics, the simulation mass density is scale- free, so we can- 
not directly read out the disc column mass density in physi- 
cal units. We dimensionalize the simulation density by solv- 
ing for E in the vertically integrated mass conservation equa- 
tion after having picked a constant accretion rate M:* 



(8) 



where r is the BL radial coordinate, and is the mass- 
averaged radial velocity for the GRMHD simulations, com- 
puted via: 



J 8=0 



(9) 



4-3.1 Radial velocity matching (a) 

At large radii (r > Tie), the simulation is no longer con- 
verged, and we resort to matching the simulation result with 
an analytical disc model. Rather than directly using the 
Novikov & Thorne (1973) disc model, we re-solve the NT 
disc equations using our matched GRMHD luminosity pro- 
file (see Appendix B for details). The final matched radial 
velocity profile is shown in Figure 5. 

We choose the disc viscosity a so as to ensure conti- 
nuity in the radial velocity profile (Fig. 5). The matching 
values of a for each disc are listed in Table 1. The loca- 
tion of the radial velocity matching point is set to be the 
same as the luminosity matching point. To verify that our 
Q value from radial velocity matching is sensible, we com- 
pare it to the directly computed value from the simula- 
tion Qsim = T.7/{hPtot) where r-7 is the height integrated 



° We set M to be the value corresponding to the disc luminosity 
as determined in §4.1.1. They are related by L = r]Mc^ where 
via*) is the spin dependent accretion efficiency of the NT disc. 
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Figure 5. A comparison between the mass averaged simulation 
radial velocity profile (dotted) and the final matched radial veloc- 
ity profile (solid line) . The matching point is depicted as the large 
coloured circle, whereas the vertical dashed line denotes the loca- 
tion of the ISCO. The gaps in the GRMHD radial velocity profile 
represent grid cells whose radial velocity is pointed outwards. 
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Figure 6. The vertically averaged a profile as computed from 
the GRMHD simulation (dotted) compared to the a obtained 
by the radial velocity matching method (solid horizontal lines, 
corresponding to the a values listed in Table 1). For reference, we 
have plotted the ISCO locations using dashed vertical lines. 



shear stress, and /iPtot is the height integrated total pres- 
sure (where Ptot ~ Pmag + Pgas). The near continuity in the 
simulation a profile with the matched value suggests that 
our method for determining a is sound. 

From the radial velocity profile and from our choice of 
constant AI , the mass continuity equation (Eq. 8) yields the 
final matched disc column density profile (Fig. 7). 



5 RAY TRACING 

The final step in the process of generating the accretion 
disc's SED is to sum the light from all the individual annuli. 



Figure 7. Dimensioned column density profiles for the GRMHD 
discs (solid lines) compared to NT discs (dotted lines). Note the 
sharp drop in NT column density as the fluid reaches the plunging 
region, in contrast to the more gradual tapering observed in the 
GRMHD simulations. 



The observed disc image is determined by shooting a series of 
parallel light rays from an image plane that is perpendicular 
to the black hole line of sight. We partition the image plane 
into a sequence of polar grid cells with logarithmic radial 
spacing (to concentrate resolution near the black hole), and 
with uniform angular spacing. To produce the ray traced 
images, we shoot rays for Nr = 300 radial grid cells and 
A^0 = 100 polar cells. The grid is then further squashed by 
a factor of cosi to match the geometry of the inclined disc. 

A single ray is shot for each grid cell, and by numerically 
integrating the (second-order) geodesic equations for light: 



where A is an affine parameter along the geodesic, and FJ^^ 
are the connection coefficients, we locate the first disc mid- 
plane intersection for each light ray (we stop at the disc 
midplane since we have an optically thick disc - see Figure 
11 for a plot of the optical depth profile). For simplicity, 
we opt to use the disc midplane as the point of light emis- 
sion rather than using each annulus' precise photosphere. 
K2011 suggested that the latter is important only for thick 
discs viewed at high inclinations angles (which is when disc 
self-shadowing becomes important). In this work, we do not 
expect the distinction between midplane and photosphere 
to be important since we are dealing with thin discs (where 
h/r <^ cosi). The ray tracing code used in this work was 
originally developed by Scherbakov & Huang (2011), with 
further refinements by K2011. 

For each ray, we compute (by interpolation on the 
GRMHD grid cells, see Appendix C for details) the local 
values of: flux T^a, column density E, vertical gravity Q, 
and comoving light ray incidence angle /i. The local comov- 
ing spectrum is then obtained by interpolating on our grid 
of TLUSTY atmosphere models. We then apply the relevant 
Doppler boosting and gravitational redshifting to these spec- 
tra (see Fig. 8 for the ray traced images) . Finally, the overall 
disc spectrum is obtained by integrating the spectra corre- 
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Figure 8. Colour images of the inner disc (r < 15M) produced via ray tracing (viewed at an inclination angle oi i = 60°) for models A, 
B, and C in Table 1. The colours correspond to the flux Fi, integrated over different energy bands, where red colour corresponds to the 
energy band E < 4keV, green 4keV < E < 12keV, and blue E > 12keV. For each value of BH spin, the colour mapping for the GRMHD 
and NT panels are identical, and each colour channel is normalized to the maximum flux in that channel. Note the appearance of a 
blue spot in the GRMHD disc plunging region. This blue spot is responsible for the appearance of a nonthermal high energy power-law 
feature in the disc spectra (see Figure 9). 



spending to these light rays over the apparent disc area in 
the image plane. We take a fiducial BH distance of 10 kpc 
when generating the disc spectra. 



would be an asymmetric error in BH spin estimates (since 
the current fitting models are based on the cooler NT disc). 
The quantitative size of this spin bias is discussed in §6.2. 



6 RESULTS 

Qualitatively, we can spot a few differences in Figure 8 be- 
tween the GRMHD (top) and NT (bottom) panels. We see 
that for spinning black holes, the plunging fluid emits in 
the hard X-rays (appearing as a > 12keV blue smudge in 
Fig. 8). Even outside the plunging region, there is a notice- 
able increase in disc luminosity, which results in harder an- 
nuli spectra everywhere (compare the bright white Doppler 
spot in the upper right (Model C, GRMHD) panel, with the 
duller orange spot in the lower right (Model C, NT) panel. 

By integrating the flux in the image plane, we arrive at 
the observed disc spectrum (see Fig. 9). For spinning black 
holes, the spectra corresponding to the GRMHD simulations 
appear to exhibit a power-law component at high energies. 
We discuss this effect in §6.1. In addition, for all models, the 
location of the energy peak in the spectra for the GRMHD 
discs is shifted towards higher energies. The harder overall 
spectrum of GRMHD compared to NT implies that there 



6.1 Power law tail 

Observationally, high-energy power-law tails are often 
seen in the X-ray spectra of black hole binary sys- 
tems (e.g. Miyamoto et al. 1991; McClintock et al. 2001; 
Reis et al. 2010). The current leading theory to ex- 
plain these high-energy power-law phenomena is the ac- 
tion of a hot optically thin corona (Zhang et al. 1997b), 
which produces the power-law tail by means of ei- 
ther Comptonization of disc photons, synchrotron radia- 
tion, and/or synchrotron self-Comptonization (Miller 2007; 
McClintock & Remillard 2006). Despite the many proposed 
coronal models (e.g. Haardt & Maraschi 1991; Dove et al. 
1997; Kawaguchi et al. 2000; Liu et al. 2002), there is no 
agreed-upon standard (i.e. there is scarce agreement on the 
geometry and physical properties of the corona). Although 
a hot corona is often invoked to explain the existence of a 
high energy power-law tail, in this work we naturally re- 
cover a hot power-law tail from just the thermal disc fluid 
that occupies the plunging region. 



Light from 



I = 60° 




10-1 jqO 

E [keV] 



Figure 9. GRMHD and NT disc spectra for eacii of the 6 mod- 
els listed in Table 1. The flux normalization corresponds to an 
assumed distance of 10 kpc. Note that the GRMHD spectra peak 
at slightly higher energies than their NT counterparts. Note also 
the emergence of a high energy power-law tail for the spinning 
black holes. The dashed spectra (labelled GRMHD2) correspond 
to the alternate GRMHD luminosity profile (dashed line in Figure 
1), discussed in §7.2 and Appendix D. 



To better understand the origin of the apparent high- 
energy power-law tail in our models, we examine the area- 
weighted local spectrum for several annuli in model B (a* = 
0.7), chosen since this simulated disc exhibits the strongest 
tail relative to the disc continuum. From Figure 10, we see 
that the high energy power-law tail stems solely from the 
emission of annuli inside the plunging region (coloured in 
blue). Furthermore, we see that the power-law emerges from 
the combined emission of plunging annuli at various loca- 
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Figure 10. Local emergent fluxes from various annuli for model 
B (a, = 0.7), weighted by their face-on emitting areas (i.e. A = 
27rrAr, where Ar is the radial width of the annulus). Plunging 
region annuli are coloured blue, whereas annuli beyond the ISCO 
arc coloured red to black (for reference, the ISCO is located at 
3.4Af, and the horizon is at 1.7M). Doppler and gravitational 
redshift corrections have not been applied to these spectra. 

tions (i.e. the envelope of spectra from r = 2M — 3M com- 
prise the power-law) . No single annulus emits the full power- 
law. 

In our plunging region model, we find that the strength 
of the high energy power-law also depends on BH spin. The 
spinless models (C, D, E, F) do not appear to produce any 
significant power-law component. The reason for why model 
B (a* = 0.7) produces a stronger power-law than model A 
(a* = 0.9) is simply that the plunging region in B covers a 
larger area in the image plane than A, leading to a larger 
overall plunging region flux (compare the relative sizes of 
the blue plunging region blobs in Figure 8). We also find a 
correlation in the strength of the power-law with disc view- 
ing angle; edge-on discs produce the strongest power laws 
relative to the thermal continuum since Doppler beaming 
increases with inclination angle, and hence acts to preferen- 
tially boost the plunging region emission. 

We find that the spectral hardening factor for the plung- 
ing annuli appears to start growing non-linearly as the gas 
approaches the horizon (see Fig. 10 and note the rapid shift 
in the spectral peak locations for the blue plunging region 
annuli, despite that T^g is nearly constant in the plung- 
ing zone as per Fig. 4). The physical reason for this rapid 
spectral hardening is simply that the plunging gas quickly 
becomes extremely hot. Upon entering the plunging region, 
the fast radial plunge causes the annuli to rapidly thin out 
in column density. The dropping column mass allows the 
annuli to make the transition from being optically thick to 
effectively optically thin^ (see Fig. 11). Despite this drop in 
column density, the gas needs to maintain a roughly constant 
cooling rate (i.e. Tcff remains constant in the plunging zone - 
see Fig. 4). A property of effectively optically thin gas is that 

^ An effectively optically thin medium is one that is optically 
thin to absorption, and optically thick to scattering. 
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Figure 11. Optical depths corresponding to model B (a* = 0.7). 
The total optical depth m = T,kji is computed using the Rosse- 
land mean opacity kji (including both absorption and scattering) . 
The effective optical depth is computed as r^g = ^/Srabsfscat 
where Tabs and Tgcat are the individual absorption and scattering 
optical depths respectively. 



it cannot radiate efficiently (Stiapiro, Ligtitman, & Eardley 
1976). To fceep a constant cooling rate despite the rapidly 
dropping optical depth, the gas must heat up tremendously, 
leading to very hot emission. 

The power laws produced by the plunging region in 
our models have photon power-law indices that range from 
r — 4 — 6, with scattering fractions^" ranging from 1-15 
per cent. In our model, we find that the scattering frac- 
tion (and hence strength of the power-law) increases mono- 
tonically with inclination angle, whereas the photon power- 
law index decreases monotonically. The scattering fraction 
range spanned by our models is comparable to the observed 
range covered by the thermal-dominant state, as well as 
much of the range of the intermediate, and steep power- 
law states for black hole binaries (which ranges from 0-25 
per cent - see Steiner et al. 2009b, 2011). However, even in 
the steep power-law state, the observed photon index range 
of r = 2.2 - 2.7(Gou et al. 2009, 2011; Steiner et al. 2011) 
is much less steep than what we obtain from the plunging 
region of our models. This discrepancy may simply be due 
to a limitation of our model, which neglects the physics of 
disc self-irradiation. Including this effect would result in a 
hotter disc that would act to boost the strength of the high 
energy tails. Also, self-irradiation becomes increasingly im- 
portant as one approaches the black hole (where light bend- 
ing is most severe). This implies that including irradiation 
will likely result in power laws that are less steep (since the 
emission from the hottest, innermost annuli will be boosted 
the most by this self-irradiation effect). We speculate that 
a more complete treatment of the plunging region (that in- 
cludes the physics of disc irradiation and magnetic fields) 



would produce high energy tails that are better matched to 
observations. 

Since many of our approximations become poor in the 
plunging region (see §7 for a discussion of some of the pit- 
falls of our method) , we strongly caution against reading too 
much into the quantitative details of the power laws of Fig- 
ure 9 (i.e. slope, normalization). The high energy power-law 
result should be taken only at the qualitative level; although 
the plunging region luminosity is dwarfed by the disc lumi- 
nosity, it starts to dominate the flux at i5 > 20 keV taking 
on a nonthermal shape that appears as a power-law. 

6.2 Quantitative effect on spin 

As already discussed, the thermal emission from the disc is 
slightly stronger and hotter in the GRMHD models com- 
pared to the equivalent NT model. This will introduce an 
error in BH spin estimates obtained from the continuum fit- 
ting method. To quantify the effect, we generate mock ob- 
servations of our GRMHD discs and we fit these simulated 
observations with the currently used suite of continuum- 
fitting disc models. For simplicity, to generate these mock 
observations, we use an idealized irrstrument with effec- 
tive area that is independent of photon energy. We choose 
the instrument's effective area'^^ to be 1000 cm^, and use a 
3h exposure time. We compute photon statistics in 1000 
energy chanrrels uniformly spaced in log E, rangirrg from 
0.4 keV < < 50.0 keV. 

We then fit our simulated disc observations with bh- 
SPEC (Davis et al. 2005; Davis & Hubeny 2006a), which is 
a set of (TLUSTY based) disc spectra often used in con- 
tinuum fitting. The bhspec spectra are computed in much 
the same way as our procedure in §1, except that bhspec 
uses a classic NT disc instead of our GRMHD disc. When- 
ever necessary, we also attempt to fit a power-law tail to the 
models using SIMPL (Steiner et al. 2009a), which is a physi- 
cally motivated Comptonization model that has been quite 
successful in fitting many observed high energy power-law 
tails (Steiner et al. 2009b). 

The fitting is handled through xspec, a spectral fit- 
tirrg software package commonly used by X-ray astronomers 
(Arnaud 1996). For each spectrum, we fit for only two pa- 
rameters: BH spin, and mass accretion rate. We fix the 
mass, distance, and inclination to exactly the values used 
to generate the simulated observations. In the cases where a 
power-law fit via SIMPL is needed, we fit for two additional 
power-law parameters (essentially the normalization and the 
spectral slope of the power- law, which SIMPL errcapsulates as 
the scattering fraction of soft disc photons and the photon 
power-law index). We list in Table 2 our spin results from 
this fitting exercise. The first three columns correspond to 
the following: 

'NT' - We perform a bhspec fit on spectra computed 
from a NT disc model, with disc parameters corresponding 
to Table 1. We use this as a baseline for comparing with 
our GRMHD fit results. For the fit, we use an energy range 



The scattering fraction denotes the fraction of thermal seed 
photons that undergo Comptonization to produce the power-law 
component. 



This is modelled after the Rossi X-Ray Timing Explorer's 
total effective area in the hard energy band (> 20keV), see 
Jahoda et al. (2006) 
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NT GRMHD GRMHD+PL GRMHD2t GRMHD2+PLt 

(Model A - a. = 0.9, L/L^dd = 0.35, a = 0.22) 



15 


0.902±0.004 


0.909±0.018 




0.913 ±0.017 




30 


0.909±0.005 


0.914±0.015 




0.916 ±0.016 




45 


0.904±0.006 


0.912±0.013 




0.916 ±0.013 


0.916±0.005 


60 


0.902±0.005 


0.919±0.012 


0.913±0.006 


0.921 ±0.011 


0.911±0.004 


75 


0.897±0.007 


0.925±0.011 


0.916±0.005 


0.934 ±0.010 


0.908±0.004 



(Model B - a, = 0.7, L/L^dd = 0.32, a = 0.10) 



15 


0.68±0.01 


0.73±0.03 


0.72±0.03 


0.75±0.03 


0.73±0.03 


30 


0.69±0.01 


0.74±0.03 


0.72±0.02 


0.75±0.03 


0.72±0.03 


45 


0.69±0.01 


0.75±0.02 


0.71±0.03 


0.76±0.03 


0.72±0.03 


60 


0.70±0.01 


0.75±0.02 


0.71±0.03 


0.76±0.02 


0.71±0.04 


75 


0.70±0.01 


0.76±0.01 


0.69±0.04 


0.77±0.03 


0.69±0.04 



(Model C - a, = 0, L/LEdd = 0.37, a = 0.04) 



15 


-0.12±0.04 


-0.04±0.09 




-0.01±0.09 




30 


-0.10±0.04 


-0.03±0.09 




0.01±0.09 




45 


-0.09±0.04 


0.01±0.08 




0.03±0.08 




60 


-0.08±0.04 


0.03±0.08 




0.07±0.09 




75 


-0.06±0.05 


0.05±0.09 




0.10±0.10 





(Model D - a, = 0, L/L^dd = 0.70, a = 0.04) 



15 


-0.19±0.02 


-0.01±0.09 




0.04±0.09 




30 


-0.17±0.02 


0.00±0.09 




0.05±0.09 




45 


-0.15±0.03 


0.01±0.08 




0.08±0.09 




60 


-0.14±0.03 


0.03±0.08 




0.11±0.09 


0.08±0.10 


75 


-0.12±0.04 


0.08±0.09 


0.04±0.09 


0.17±0.09 


0.07±0.13 



(Model E - a* = 0, L/Ljsdd = 0.71, a = 0.08) 



15 


-0.08±0.03 


0.04±0.08 




0.07±0.08 




30 


-0.08±0.03 


0.05±0.08 




0.10±0.08 




45 


-0.07±0.02 


0.09±0.07 




0.12±0.09 




60 


-0.06±0.02 


0.10±0.08 




0.16±0.08 




75 


-0.06±0.03 


0.13±0.09 




0.20±0.09 


0.11±0.07 



(Model F - a, = 0, L/L^dd = 0.36, a = 0.03, 1-loop) 



15 


-0.13±0.07 


-0.07±0.14 


30 


-0.11±0.07 


-0.05±0.13 


45 


-0.10±0.06 


-0.02±0.13 


60 


-0.09±0.06 


-0.01±0.13 


75 


-0.09±0.06 


0.01±0.13 



-0.03±0.15 
-0.01±0.14 
0.02±0.14 
0.06±0.13 
0.10±0.14 



0.09±0.13 



Note - All spectra were fit with BHSPEC disc spectra, where agt = 0.1. The quoted uncer- 
tainties correspond to the systematic GRMHD luminosity profile uncertainties (determined 
empirically by analysing the last 5 subsequent lOOOM time chunks of the GRMHD simu- 
lations). The spectral fitting statistical uncertainties were negligibly small (Aa, ~ 0.001). 
This is because the disc parameters that are responsible for most of the uncertainty (mass, 
distance, inclination) were not allowed to vary during the fit. 



t These two rightmost columns correspond to a disc model that uses an alternative measure 
of the GRMHD luminosity (See §7.2 and Appendix D for details). The luminosity profiles 
corresponding to GRMHD2 are shown in Figure 1. 
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Figure 12. NT disc spectra for model C (a* = 0) computed for 
two ciioices of a. The slight shift in the location of the spectral 
peak causes a systematic error in the recovered BH spin if Qfjt 
does not match the instrinsic a of the disc. 

of 0.4-8.0 keV. 

'GRMHD' - We perform a bhspec fit on the GRMHD 
disc spectra on the energy range 0.4-8.0 keV. We employ 
an 8.0keV upper energy cutofT in the fitted spectrum since 
we want to exclude the power-law spectral feature. 

'GRMHD+PL' - In addition to continuum fitting, we 
also fit for the power-law. Formally, we fit a bhspeC(8)SIMPL 
model over the energy range 0.4-50.0 keV (the full range of 
our mock observation). 

The primary purpose of the 'NT' column is to disen- 
tangle another systematic effect that is independent of the 
extra luminosity of the GRMHD simulations. The bhspec 
spectral model is tabulated only for am = 0.1, however in 
our disc models, only the a, = 0.7 disc happens to have a 
matching a = 0.1 (see Model B in Table 1). This mismatch 
between the fitting model ant and the input model a leads to 
a systematic bias in the recovered spin. The spectral depen- 
dence on a arises from the fact that column density scales 
inversely with a, so discs with higher a values tend towards 
lower optical depths. A lower optical depth medium has gas 
that is hotter, which results in harder spectra. In Figure 12, 
we illustrate the infiuence that the choice of a has on the 
spectum of model C (a* = 0). As expected, we find that the 
higher a disc has a harder spectrum (Fig. 12). 

In the context of spin fitting, if Qflt > ct, we expect the 
fitted spins to be too low (since for any given spin, the fitting 
spectrum will be harder than the intrinsic disc spectrum), 
and vice versa for ant < a. This effect is seen in the 'NT' 
column of Table 2 (which differs from the fitting model only 
by its choice of a). Models C, D, E, and F all have agt > a, 
and consequently the recovered spins are too low. On the 
other hand, model A has ant < a, yielding recovered spins 
that are too high. 

We interpret the difference between the 'NT' column 
and the 'GRMHD' column as the effect of the extra lumi- 



nosity from the inner regions of the GRMHD discs, which is 
to systematically increase the recovered spin. The spin de- 
viation also grows with inclination angle (see the 'GRMHD' 
column in Table 2), since the net effect of Doppler beaming 
(which becomes more important at higher inclinations) is to 
enhance the inner disc emission. We also find that the spin 
deviation becomes smaller at high spins (this is due to the 
fact that the ISCO-spin relation becomes steeper at high 
spins) . 

In those cases where using the SIMPL model significantly 
improves the fit (i.e. when the fit improved by more 
than a factor of 2, upon inclusion of a power-law fit com- 
ponent), the modelling of the power-law also leads to bet- 
ter agreement between fitted and input BH spin (compare 
'GRMHD' and 'GRMHD-fPL' columns in Table 2). By fit- 
ting a power-law component, some of the hot photons get 
associated with the power-law, which lowers the fitted spin 
(since the remaining disc spectrum looks softer after removal 
of the hard power-law photons). However, we note that the 
GRMHD luminosity excess extends even beyond the ISCO 
(see Fig. 1). This explains why even fitting for the power- 
law does not completely eliminate the upwards spin bias 
(i.e. the 'GRMHD-fPL' columns still yield spins that are 
higher than the 'NT' column in Table 2); although some of 
the excess luminosity is bound up in the nonthermal plung- 
ing region (which can be excluded by means of power-law 
fitting), there is still some residual thermal disc luminosity 
excess, that acts to bias the GRMHD spin upwards. Over- 
all, the effect of the extra plunging region light results in 
Aa* ~ 0.1, 0.06, 0.03 corresponding to spins a* = 0, 0.7, 0.9 
respectively (estimated by comparing the 'NT' column with 
the 'GRMHD' column in Table 2 for the ~30 per cent Ed- 
dington discs). 

These systematic spin errors should be thought of as 
an upper bound since the discs that we consider correspond 
to fairly high accretion rates (with L/I/Edd > 0.3). Previ- 
ous work (Paczyriski 2000; Shafee et al. 2008a) suggests that 
thinner discs (corresponding to lower accretion rates) better 
match the NT model. Our results in Table 2 also support 
this claim; in comparing the spin fits from model C (thin 
disc) with models D and E (thick discs), we find that model 
C agrees best to the NT disc (i.e. that model C has the small- 
est difference between the 'NT' and 'GRMHD' columns). In 
addition, model C also exhibits the least deviation from a 
thermal spectrum at high energies (compare the > 20keV 
emission from model C with those of D and E in Fig. 9). 

We did not simulate lower accretion rate discs due to the 
rapidly increasing computational cost associated with such 
simulations. Typically, continuum fitting is not applied to 
systems where the accretion rate exceeds 30 per cent Edding- 
ton since beyond this critical point, the continuum fitting 
method is no longer robust to variations in the disc accre- 
tion rate (see McClintock et al. 2006; Steiner et al. 2009b). 
One possible explanation is that disc self-shadowing be- 
comes important beyond this critical accretion rate thresh- 
old (Li et al. 2010). 

To put into perspective how significant our quoted spin 
deviations are, we compare our results to the spin uncertain- 
ties from actual continuum fitting exercises in the literature 
(See Table 3). Since the current observational uncertainties 
are significantly larger than the deviations that we find in 
this work, we conclude that the extra light from the inner 
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and plunging regions of the disc do not limit our ability to 
make accurate spin estimates through the continuum fitting 
method. The current limiting factors in measuring BH spin 
are how accurately one can determine the distance, mass, 
and inclination of black hole binary systems. 



7 DISCUSSION 

The disc model that we adopt in this work is not completely 
self-consistent. The annuli that we compute using tlusty 
are assumed to have reached their equilibrium configuration. 
However, close to the black hole there is simply insufficient 
time for the fluid to reach equilibrium (see Fig. 13 for a 
plot of the various time-scales). In addition, the GRMHD 
simulations show that the pressure in the innermost regions 
of the disc is magnetically dominated, yet in the tlusty 
annuli calculations, we completely ignore magnetic pressure 
support. Previous work by Davis et al. (2009) suggests that 
including magnetic support may lead to a slight hardening 
of the annuli spectrum. Another problem for the innermost 
annuli is that they become extremely hot (Tgas > 10* — lO' 
K), which may invalidate tlusty's handling of Comptoniza- 
tion (i.e. an angle averaged Kompaneets treatment, which 
is only valid for nonrelativistic electrons, see Hubeny et al. 
2001). Finally, due to the increasingly strong light bending 
effects near the black hole, the process of disc self-irradiation 
(i.e. where one part of the disc shines on another part) may 
also become important, especially for the innermost annuli 
(Li et al. 2005). The net effect of this disc self-irradiation 
would be to further heat up the annuli, causing additional 
spectral hardening. We have ignored disc self-irradiaton in 
our spectral model since this reprocessing of light cannot be 
easily handled by our grid of precomputed annuli. 

These issues primarily affect annuli that are very close 
to the black hole, which radiate in the hard X-rays (> 
20 keV) . We do not expect these inconsistent annuli to affect 
our spin fitting results in Table 2 since the recovered spin 
is mainly constrained by photons near the peak of the spec- 
trum (located at a few keV in Figure 9). The inconsistent 
annuli will however affect the high energy part of the spec- 
trum, which is why we caution against reading too deeply 
into the quantitative details of the power laws. 



7.1 Other signatures of the plunging region 

Observations of disc variability could potentially be used to 
discriminate between our plunging region model and other 
coronal models. Assuming that the variability is caused by 
Doppler beaming of hotspots moving about in the disc fiow, 
we associate the hotspot orbital time-scale with the vari- 
ability time-scale. Since the orbital time decreases mono- 
tonically as one approaches the horizon, the variability from 
the inner plunging region of the disc ought to be more rapid 
than that occurring farther out in disc. Hence, our plung- 
ing region model predicts that the hard X-ray light (origi- 
nating from the innermost regions of the disc) would have 
faster variability than the soft X-ray light (produced farther 
out in the thermal disc). Furthermore, in the power-density 
spectrum in the hard X-ray bands, we expect most of the 
variability power to occur above the ISCO orbital frequency 



Model B - (a* = 0.7) 




r/M 

Figure 13. Several time-scales of interest are plotted here for 
model B (a, = 0.7). We show the accretion time-scale (solid 
blue), the photon diffusion time-scale (dash-dotted green), the 
photon absorption time-scale (dotted red), and the photon scat- 
tering time-scale (dashed dark cyan). The absorption and scat- 
tering time-scales are computed via t = l/(pKc), using midplane 
gas densities (p) and opacities (k). The diffusion time-scale is 
computed via = T^^^^tsca.t where Tacat is the scattering optical 
depth to the disc midplane (i.e. we expect t^^^^ scatterings to oc- 
cur before a photon can diffuse out). Note that close to the black 
hole, the accretion time-scale becomes shorter than the photon 
absorption time-scale, which indicates that the fluid close to the 
black hole has insufficient time to reach thermal equilibrium. 



(since we associate this hard emission with the plunging re- 
gion). 

Polarization observations could also serve as a means 
to distinguish between other coronal models and a plunging 
region model. Generally, one expects evolution of the po- 
larized fraction and net polarization angle as the observed 
photon energy is varied. This is due to higher energy photons 
(which are typically emitted closer to the black hole) feel- 
ing stronger relativistic effects, leading to a stronger shift 
in their polarization vectors (Schnittman & Krolik 2009, 
2010). The planar geometry of the plunging region may lead 
to a polarization signature that differs significantly from 
coronal models. If the high-energy power law originates from 
the plunging region, we expect that its polarization proper- 
ties will vary continuously from the thermal component to 
the non-thermal power law. This contrasts with some alter- 
native coronal geometries (e.g. spherical models) that pro- 
vide abrupt transitions in the polarization properties from 
one component to the next (Schnittman & Krolik 2010). 

In the plunging region model, even at energies where 
the power law dominates, the variation of the polarization 
fraction and angle with energy could be consistent with 
thermal state polarization models that assume a thin disc 
geometry. The signal might then be similar to the models 
of Li, Narayan, & McClintock (2009); Schnittman & Krolik 
(2009), who calculate the polarization due to direct emis- 
sion from a NT model, assuming a (modified) blackbody 
spectrum and intrinsic disc polarization appropriate for a 
scattering dominated atmosphere. The models also include 
the effect of scattered emission from light that is strongly 
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Table 3. Spin measurements yielded by continuum fitting. 



Black hole 


1* 


in 




Qflt 


Reference 


A0620-00 


0.12 ± 0.19 


51.0 ±0.9 


0.11 


0.01+0.1'' 


Gou et al. (2010) 


H1743-322'= 


0.2± 0.3^^ 


75 ± 3 


0.03-0.3 


O.Ol+O.l'' 


Steiner et al. (2012) 


LMC X-3 


< 0.3 


67 ±2 


0.1-0.7 


0.01, 0.1*" 


Davis et al. (2006b) 


XTE J1550-564 


34+°-^ 


74.7 ± 3.8 


0.05-0.30 


0.1 


Steiner et al. (2011) 


GRO J1655-40'* 


0.70 ± 0.1 


70.2 ± 1.2 


0.04-0.1 


0.1 


Shafee et al. (2006) 


4U 1543-47'^ 


0.80 ± 0.1 


20.7 ± 1.5 


0.06-0.1 


0.1 


Shafee et al. (2006) 


M33 X-7 


0.84 ± 0.05 


74.6 ± 1.0 


0.07-0.11 


0.01 


Liu et al. (2008, 2010) 


LMC X-1 


92+0 05 


36.4 ± 2.0 


0.15-0.17 


0.01+0. 1=* 


Gou et al. (2009) 


GRS 1915+105'* 


> 0.95 


61.5-68.6 


0.2-0.3 


0.01, 0.1*" 


McClintock et al. (2006) 


Cygnus X-1 


> 0.95 


27.1 ±0.8 


0.018-0.026 


O.Ol+O.l'' 


Gou et al. (2011) 



Note — The spin uncertainties correspond to the 68 per cent (Icr) level of confidence, whereas the 
inequalities are to the 3a" level. 



^ The spin errors are marginalized over both choices of a 
^ The spin constraint covers both choices of a 
No reliable mass estimate is available for this source 

The quoted spin errors have not been rigorously computed, and have been arbitrarily doubled from the 
published estimates since these are among the first systems for which continuum fitting was applied). 



bent by the relativistic spacetime of the black hole and rein- 
tersects the disc surface (i.e. returning radiation). In these 
calculations, matter in the plunging region is allowed to scat- 
ter returning radiation, but does not provide direct emission. 
Since the intrinsic emission and returning radiation have 
different polarization signals, a similar calculation that in- 
cludes the direct (power law) emission from the plunging 
region is required to confirm the above speculation and dis- 
criminate between coronal models. The impending launch 
of NASAs Gravity and Extreme Magnetism Small Explorer 
(Swank et al. 2010; GEMS) mission promises observational 
constraints that may be capable of discriminating between 
such models. 



7.2 How does our choice of cooling function 
influence the results? 

A criticism that might be levied against the GRMHD sim- 
ulations is that our choice of cooling function (as defined in 
Eq. 1) is arbitrary. We have interpreted this cooling func- 
tion as the rate of radiative cooling, despite the fact that 
the GRMHD simulations do not inherently model radiation 
physics. 

We would like to test how our choice of cooling func- 
tion impacts the results of §6. Rather than rerunning the 
GRMHD simulations with another prescription for cooling 
(i.e. by changing the functional form of Eq. 1), we use the 
simulation's dissipative heating profile as another means to 
get a cooling rate. Instead of simply assuming that dissipa- 
tion equals cooling locally, we consider the physics of energy 
advection (i.e. that the fluid releases its heat downstream). 
To generate a disc model that is more self-consistent with the 
TLUSTY annuh, we use the TLUSTY vertical structure 
when computing the rate of energy advection (see Appendix 
D for a thorough exposition on this method of computing 
the new GRMHD luminosity profile). The new dissipation- 
based cooling function (which we label as GRMHD2) re- 



leases slightly more luminosity inside the plunging region 
(compare the solid and dashed lines in Figure 1). 

We repeat the same exercise as in §6 with this new 
GRMHD2 derived luminosity profile. The spectra corre- 
sponding to GRMHD2 are shown as the dashed lines in Fig- 
ure 9. Without fitting for the power laws, we find that the re- 
covered spins are slightly higher owing to the extra plunging 
region luminosity (compare the 'GRMHD' and 'GRMHD2' 
columns in Table 2). However, the more realistic fit (i.e. 
including the power-law component) takes care of this ex- 
tra light, yielding comparable spin estimates to before. The 
enhanced power laws are purely due to the extra plunging 
region luminosity (compare solid and dashed lines in Figures 
1 and 9). The worst cases for the 30 per cent Eddington discs 
now have Aa* ~0.15, 0.07, 0.03, for a, = 0, 0.7, 0.9 respec- 
tively. Given these results, we still conclude that the extra 
luminosity from the GRMHD inner disc does not limit our 
ability to measure BH spin through the continuum fitting 
method (the dominant source of uncertainty is still the ob- 
servational uncertainties in distance, mass, and inclination). 

As in K2011, we estimate that Aa* ^ 0.15, which is 
lower than the estimate of Noble et al. (2011), who in- 
fer Aa ~ 0.2 — 0.3 for a Schwarzschild black hole. Due to 
the reduced sensitivity of the continuum fitting method at 
low spins, this amounts to a somewhat modest discrepancy, 
which we attribute to various differences in the setup and 
intitial conditions of the GRMHD simulations used by the 
different groups. 



7.3 What is the impact of the initial magnetic 
field topology? 

It can be argued that our choice of four weak poloidal mag- 
netic loops is arbitrary, and that the most natural choice is 
either a singly looped purely toroidal or purely poloidal field 
(Igumenshchev, Narayan, & Abramowicz 2003). Our choice 
is motivated by the belief that nature does not begin the ac- 
cretion process with an ordered field spanning large spatial 
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scales; starting from a weak and disordered field, we believe 
the physics of the magneto-rotational instability (MRI) nat- 
urally selects the correct magnetic configuration, forgetting 
about the initial conditions. Ideally, we would like to start 
the disc with an infinite number of loops to model this tan- 
gled geometry, however owing to resolution limitations we 
settle with four loops. 

To determine how much influence the initial field geom- 
etry could have on disc spectra, we compare the results of 
Model C (4-loop) and Model F (1-loop). In general, we find 
that the I-loop case releases significantly more luminosity 
inside the plunging region. The large-scale magnetic fluxes 
from the I-loop geometry produce more stress and dissi- 
pation than its 4-loop counterpart. However, most of this 
energy is released at large polar angles (i.e. far from the disc 
midplane where the gas is unbound), and for the purposes 
of modelling the thermal disc emission, we choose to ignore 
this unbound gas contribution to the luminosity. 

Making this distinction between bound and unbound 
gas, we find that the luminosity profile corresponding to the 
bound gas in the I-loop case closely matches that of the 
4-loop case. This leads to very similar disc spectra (com- 
pare the spin fit results from Models C and F in Table 
2 and spectra in Figure 9). We conclude that in the con- 
text of the thermal continuum for thin discs, the resul- 
tant spectrum is largely independent of the initial mag- 
netic field geometry. This is in contrast to the case of 
magnetically arrested accretion flows, whose final accretion 
state and jet power depends crucially on both the strength 
and geometry of its seed magnetic field (Igumenshchev 
2008, 2009; Tchekhovskoy, Narayan, & McKinney 201 1; 
McKinney, Tchekhovskoy, & Blandford 2012). For thin- 
discs, the primary differences arising from the choice of I- 
loop and 4-loop configurations are highlighted in §8 and 
§9 of Penna et al. (2010). For a complete discussion of 
this topic, see also: K20II; Noble et al. (2010, 20II); and 
Hawley, Guan, & Krohk (20II). 

7.4 Convergence of GRMHD simulations 

Recent work by Hawley, Guan, & Krolik (20II) suggests 
that the growth and saturation of the MRI in disc simu- 
lations is highly dependent on their vertical and azimuthal 
resolution. These authors propose the following two conver- 
gence criteria for the MRI to be well resolved: the vertical 
criterion is given by Qz — Amri/(A2:) > 10, and the az- 
imuthal criterion Qy = Amri/(Aj/) > 5, where Amri is the 
MRI length scale in the direction that Q is being measured. 

To compare with Hawley, Guan, & Krolik (20II), we 
adopt their prescription for Amri = 2-kH l3~^^'^\Bz\/\B\ in 

— 1/2 

the evaluation of Qz and Amri = 2-KHj3y in Qy, where H 
is the disc scale height, plasma P = -Pgas/-Pmag, and plasma 
Py is computed with only the toroidal field component. For 
the 4-loop standard resolution runs^'^, we find = 6 — 13 
at the disc midplane and Qy = 2 — 3. 

As an alternative model-free means to estimate the 
MRI length scale, we have also performed a correlation 
length analysis similar to the one presented in §3.8 of 



McKinney, Tchekhovskoy, & Blandford (2012). From the 
spatial autocorrelation function for gas density, we find 
Qz ~ 9 — II and Qy ~ 5 — 7. Azimuthal shear elongation 
of the MRI unstable zone is the reason why the correlation 
length based Qy is larger than the Hawley, Guan, & Krolik 
(20II) estimate for Qy. 

Since our value for Qy is low, this implies that the sim- 
ulations are only marginally resolving the azimuthal disc 
substructure (and hence MRI turbulence). A consequence 
of this low resolution is that the strength of the MRI turbu- 
lence (as measured by the disc viscosity a) does not appear 
to converge with respect to resolution. Comparing models 
D and E in Table I, we find that a has doubled upon dou- 
bling both cj> and 6 resolutions. However, recent work by 
Sorathia, et al. (2012) suggests that q is a poor metric for 
gauging simulation convergence since large variations in a 
axe observed even in the cases of abundant resolution. 

For the purposes of spectral modelling we ultimately 
only care about the simulation luminosity profile. Despite 
the low azimuthal resolution of the Penna et al. (2010) sim- 
ulations, the luminosity profile remains fairly insensitive to 
changes in ^-resolution. The convergence tests performed by 
Penna et al. (2010) show that the luminosity profile changes 
only by 10%, even after a 4- fold increase in (^-resolution. 
This robustness in luminosity proflle is also borne out when 
comparing the spectra of Model D and Model E (see Figure 
9 - these two runs only differ in grid resolution). 



7.5 Equation of state 

Another inconsistency in the GRMHD simulations is that 
we adopt an ultrarelativistic gas equation of state (with 
r = 4/3), motivated by the idea that the disc energy bud- 
get is dominated by photons (a F = 4/3 fluid). However the 
TLUSTY annuli calculations suggest that the combined pho- 
ton and gas entity acts more like a F = 5/3 fluid, especially 
in the gas pressure dominated plunging region. Despite the 
difference in F, Noble et al. (2009, 2010, 20II) have also per- 
formed GRMHD simulations of thin discs with F = 5/3, and 
they find similar results to the F = 4/3 discs of Shafee et al. 
(2008b); Penna et al. (2010). 

In figure 14, we examine the relative impact that this 
equation of state inconsistency has on the disc's entropy 
profile. A fully self-consistent model would have identical 
GRMHD and tlusty gas entropy profiles. It appears 
that well outside the ISCO, the assumption that the gas 
cools towards a fixed entropy is supported by the radiative 
transfer calculations of tlusty (Compare the fiat plateaus 
in Figure 14). However Figure 14 shows that the GRMHD 
(solid) and TLUSTY (dotted) discs reach different constant 
entropies. The alternate disc cooling model (dashed) dis- 
cussed in §7.2 yields annuli which are more self-consistent 
with the simulation entropy profile (compare dashed lines 
and solid lines in Figure 14, which now track each other 
fairly well outside the ISCO). 



This includes models A-D in Tabic 1. Model E, which has 
double the resolution, will have double the Qz and Qy values. 
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Figure 14. GRMHD vs. tluSTY vertically averaged specific gas 
entropy profiles (relative to the entropy at the horizon Shor)- The 
dimensionless entropy is computed using an ideal gas equation 
of state, where s = (^^^ ■ ln(pgas/p'"). The TLUSTY lines 
(dotted) denote the gas entropy as computed from the TLUSTY 
annuli vertical structure. The TLUSTY2 lines (dashed) represent 
annuli that comprise the GRMHD2 disc model (described in §7.2 
and Appendix D). 



8 SUMMARY 

The primary goal of this work is to determine how important 
the neglected light from the plunging region is in the context 
of BH spin measurements. To answer this question, we rely 
on GRMHD disc simulations of Penna et al. (2010), which 
capture all but the radiation physics for magnetized flow 
around a black hole. We seek to convert these dimensionless 
simulations into a form that can be directly compared with 
observations, namely the accretion discs' X-ray continuum 
spectra. 

To do this, we apply radiative transfer post-processing 
on the simulated discs. This is an advance over previous work 
on computing GRMHD disc spectra (i.e. K2011), which as- 
sumes (modified) blackbody annuli spectra. We slice the 
GRMHD disc into many individual annuli, and for each an- 
nulus, we apply a 1-dimensional radiative transfer calcula- 
tion to solve for its emergent spectrum. We sum up this col- 
lection of local spectra by means of ray tracing to get the full 
disc spectrum, which we then compare with contemporary 
disc models used to measure black hole spin. The following 
are the key results from this work: 

(1) The GRMHD based accretion discs have hotter 
spectra than the standard NT discs. The GRMHD discs 
produce more luminosity everywhere, and the contrast 
becomes most apparent inside the ISCO (see Figure 1). 

(2) The increased luminosity of the GRMHD discs com- 
pared to the classic NT discs induces a modest systematic 
bias in the derived spins of these GRMHD discs. For 
black holes of spin a, = 0, 0.7, 0.9 the spin deviation is 
Aa« ~ 0.15, 0.07, 0.03 in the worst cases (corresponding 
to inclination angles of 75°). We remark that these errors 
are well within observational uncertainties (i.e. from not 



precisely knowing the system's mass, inclination, and 
distance - see Table 3). 

(3) Without needing to invoke an external corona, the 
GRMHD discs around spinning black holes exhibit a weak 
high-energy power-law tail (Fig. 9). This power- law tail 
arises from the combined emission of the hot plunging 
region gas. The strength of this plunging region power-law 
increases with the system's inclination angle. 

In our spectral modelling approach, we have made many 
simplifications and assumptions (see §7), some of which may 
be incorrect close to the black hole horizon. Due to these 
problems, we trust the power laws presented here only to 
a qualitative level - we only trust that the plunging region 
fluid emits at much higher energies than the disc, which 
adds a nonthermal component to the overall spectrum at 
high energies. 
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APPENDIX A: LUMINOSITY MATCHING 
MODEL 

The functional form of the matching flux profile is given be- 
low (see Kulkarni et al. 2011, Page & Thorne 1974 - here- 
after PT): 

fir) = fpTir) - -J-pTf nrtN9/r^ , /pT(ne) 



[(Et _ nLt)2/n,^]^ 



C, 



when r > Tie (Al) 



_(Et - nLt)2_ 

where /(r) = 4TvrFcoin(r)/ M is the dimensionless flux, L^, 
are the specific angular momentum, energy-at-infinity, 
angular velocity (given by Eqs. 15f-h of PT), /pt is the flux 
from the Page & Thorne model (Eq. 15n of PT), and C is a 
free parameter that is determined by the torque at the inner 
disc boundary. 

The task now is to flnd an appropriate value of C such 
that the GRMHD flux and Eq. Al match at r^e. For simplic- 
ity, we factor out the radial dependence from the last two 
terms of Eq. Al to get: 



/M = /pT(r) 



(El - nvy 



(A2) 



where K = {[{El - QL"^)^ /n,r]riJpT{rie) - C} is just an- 
other constant. To ensure continuity in the flux proflle at 
Tie, we require K to be: 

■(fit _ni,t)2- 



K = 



[/GRMHD(rie) — ./'pT(»-ie)] . (A3) 



Thus, our final matched fiux profile is given by 
1 dLsLir) 



crsBT'cfr(r) 



—Aivrut 
M 
4:Tvr 



dr 



fir), 



r ^ne, (A4a) 
r > r-,e, (A4b) 



where I/bl('") is the Boyer-Lindquist luminosity profile mea- 
sured from the GRMHD simulations, and /(r) is the dimen- 
sionless flux as shown in Eq. A2. 



APPENDIX B: GENERALIZED NOVIKOV & 
THORNE MODEL 

To get the generalized NT column density and radial veloc- 
ity profile far out in the disc (to extend the GRMHD simu- 
lated disc beyond ru), we solve the following set of vertical 
structure equations (with a one-zone model for the vertical 
structure): 

(i) Vertical Pressure Balance: 

dPtot 



dz 



P9±, 



has the vertically integrated form of: 

= pQh, 



(Bla) 



where h is the vertical scale height, and Q = g±/z is the 
prescription for the local vertical gravity (as defined in Eq. 
7). 



(ii) Radiative Transfer: 

From the second moment of radiative transfer equation 
(see §1.8 of Rybicki & Lightman 1979): 



dr 



frad 
C 



we transform the dr optical depth differential to a col- 
umn mass differential dm via dr = k • dm, where k is 
the Rosseland mean opacity. We integrate the column mass 
from the surface to the disc midplane (from m = to 
m = mtot = S/2). We also assume that the flux profile 
linearly increases with mass away from the midplane (i.e. 
F{m) = Ftot (1 — m/mtot) ~ this linearity assumption is also 
used in tlusty). The vertically integrated radiative diffu- 
sion equation thus becomes 



aJ2_ _ FtotKT, 
3 ^ 4c ' 



(Bib) 



where a is the radiation constant, and T is radiation 
temperature. Note: The factor 2 discrepancy in Eq. Bib 
compared to the standard prescription for radiative diffu- 
sion stems from the linearity assumption in F(m). 

(iii) Stress: 

We adopt an oi prescription for the stress where t^- = 
Qptot. Vertical integration yields 



W ; 



ij,*dz = 2/iaPtc 



(Blc) 



(iv) Viscous Heating: 

Through the energy equation for viscous heating, it is pos- 
sible to link the heating fiux with the vertically integrated 
stress (c.f. 5.6.7-12 of Novikov & Thorne 1973). The result- 
ing expression is 



Ftot = -9.kRF{r)W, 



(Bid) 



with the dimensionless relativistic factor Rpij-) defined as 



Rpir) 



1 - 2/r, + al/rl 
1-3/r* +2a,/rl 



:V2 ■ 



(v) Equation of State: 

We ignore magnetic pressure in this analysis, and adopt 
an ideal gas equation of state, yielding 



_ pkBT aT^ 

-ftot — 1 —, 

fi 3 

where /i is the mean particle weight of the fiuid. 



(Ble) 



(vi) Column Density 

In this one-zone model, the column density represents the 
quantity 



pdz — 2hp. 



(Blf) 



Our goal is to solve Eqs. Bl with unknowns 
(Ptot,p, r, S, ft), given the values of {k, p,a, Ftot,Q, M,r). 
After some algebra, we obtain an expression for E that solves 
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Eqs. Bl. We find E by solving for the real root of the fol- 
lowing polynomial in x: 

1/4 



Model B - (a, = 0.7) 



Fi 



F1F2 



1/2 



(F2) 



10 n 
X =0, 



(B2) 



where x = E^''* and the polynomial coefficients F\,F2,Fz 
depend only on the given values: 



W 



Fi = — 



4R. 



F2 = 



2^ 



1/4 



^2p2 ■ 
^ ^tot 



(B3a) 



(B3b) 



(B3c) 



We then solve Eq. B2 at different radii to obtain E(r) 
for (r > Tie). We pick k — Kea — 0.3383 cm^ g~^ and ^ = 
0.615 tuh, which corresponds to a fluid composition of 70 per 
cent Hydrogen, 30 per cent Helium by mass. Given E(r), we 
get the mass averaged radial velocity by solving for vJ' in 
the mass conservation equation M = 2nrT,u^ . 



APPENDIX C: INTERPOLATION METHODS 

All interpolated quantities are computed by linear interpo- 
lation in log space (i.e. for two variables x, and y, we in- 
terpolate linearly on log(2;) and log(y)). We have chosen 
this interpolation scheme since it can perfectly capture all 
power-law scalings. For quantities obtained by interpolating 
on TLUSTY annuli (such as the annuli vertical structure), we 
employ trilinear interpolation over the 3 annuli parameters 
log(rcfi), log(E), and log(Q) (i.e. each parameter is linearly 
interpolated seperately in log space). 

The interpolation of spectra is handled by a more 
complicated method to account for the shape of the Planck 
function. For an optically thick medium with scattering, the 
resultant spectrum takes on the form of a modified black- 
body (c.f. Eq. 2). Since the spectra doesn't depend very 
sensitively on E, or Q, we apply the log-linear interpolation 
discussed above for these two annuli parameters. For Tctt, we 
use a more complicated interpolation scheme. Rather than 
follow Davis et al. (2005) (they applied a linear interpola- 
tion on the brightness temperature, computed with a fixed 
spectral hardening factor of / = 2.0), we switch between 
three different interpolation methods that each work for all 
choices of /. The three methods are applied to different fre- 
quency ranges of the spectrum. 

• Method 1: At low frequencies, we have that oc Tcff 
in the Raleigh- Jeans tail, and thus we can apply the linear 
interpolation method. 

• Method 2: At high frequencies, we apply linear inter- 
polation on 1/Teff and log(J^). The motivation is that for a 
blackbody-like spectrum in the high-frequency Wien limit, 
the specific intensity as a function of temperature scales as 

log(/„) oc -1/Teff. 




10 



10° 10^ 

E [keV] 



10^ 



Figure CI. Comparison of annuli spectra obtained through the 
interpolation method (solid) and an exact calculation (dashed) for 
the spin a» = 0.7 disc. The black spectra represent annuli outside 
the plunging region, whereas the red spectrum corresponds to a 
plunging region annulus. 



• Method 3: For intermediate frequencies, we take a non- 
linear combination of the two interpolation results so that 
the transition between interpolation methods is smooth. 

Denoting the interpolated intensities using methods 1 and 2 
as Iiiv) and hii') respectively, the expression for the final 
combined interpolation method is: 



exp 



ln{u2) — ln(!/) 
ln(!/2) - ln(!/i) 



In [hi u)] 



ln(z^2) - In(i^i) 



ln[l2{u)] 



Vl < I' < V2, 
V V2- 



We have set v\ = 0.07t'max, and V2 = 2.3i'max, where 
i^max is the frequency that maximizes in the spectra 
of the lowest Tcff annuli in the bracketing set used for 
the interpolation. These frequency boundary values were 
determined by minimizing the interpolation error for a pure 
Planck function. For a perfect modified Planck function 
with constant /, the chosen v\ and V2 correspond to a 
maximum interpolation error of 0.08 per cent in Ii, for 
our grid with temperature spacings AlogjQ(Tcff) = 0.1. In 
practice, the interpolation error after interpolating over the 
3 annuli parameters and emission angle is ~1 per cent (see 
Fig. CI), except for annuli in the plunging region, where 
the flux error grows to ~10 per cent (the assumption of a 
constant / is violated for annuli in the plunging region). 
The main advantage of this interpolation method is it 
provides interpolated spectra that peak very close to the 
correct frequencies, and the method works equally well 
across all spectral hardening factors. 
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APPENDIX D: AN ALTERNATIVE GRMHD 
LUMINOSITY PROFILE 

We wish to generate a disc luminosity profile that has the 
following desirable properties: 1) it is based on the GRMHD 
dissipation profile rather than cooling profile, and 2) it pro- 
duces a disc whose entropy profile is self-consistent with the 
TLUSTY annuli. We accomplish this task by considering 
the physics of energy advection. The idea is to solve for the 
radiative cooling rate, taking into account both the dissipa- 
tive heating rate (measured from the GRMHD simulations) 
and the heat advection rate (computed using the vertical 
structure from TLUSTY annuli). The energy balance equa- 
tion relating these three quantities is: 



9adv — 'Jhcat 



(Dl) 



where (Jadv is the advected heating rate per unit volume, 
9heat the viscous heating rate per unit volume, and Qcooi 
is the cooling rate per unit volume (assumed to be purely 
radiative - i.e. we ignore conduction and convection). We 
take the convention that all hatted quantities are measured 
in the comoving frame of the fluid. The advective heating 
rate is obtained by evaluating ijadv = pT{ds/dT), where p 
is the rest mass density, T is the gas temperature, s is the 
specific gas entropy of the fluid, and r is the fluid's proper 
time. Making the approximation that we have stationary 
axisymmetric flow that is devoid of motion perpendicular to 
the midplane, we have ds/dr = u^ds/dr via the chain-rule 
{ds/dr and are evaluated in BL-coordinates). Vertically 
integrating Eq. Dl while making use of (jadv ~ pTu^ (ds / dr) 
yields: 



pTu^ 



dz. 



(D2) 
(D3) 



where Qadv is the vertically integrated heat advection rate 
(evaluated with TLUSTY vertical structure), Qheat is the 
vertically integrated viscous heating rate (measured from 
the GRMHD simulated discs), and Qcooi is the net radiative 
cooling flux that escapes. For simplicity in Eq. D3, we adopt 
the following constant mass-averaged radial velocity: 



The 'sim' superscript is used to denote quantities derived 
solely from the GRMHD simulations (i.e. these quantities 
are independent of TLUSTY). Analogous to Eq. D3, the 
vertically integrated^^ GRMHD advective heating rate is 
obtained by 



ds 



(D6) 



For simplicity, we first apply azimuthal and time averag- 
ing to all simulation quantities used in Eq. D6. Since the 
GRMHD simulations are dimensionless {k/ p = 1) and em- 
ploy an ideal gas equation of state, we have that T = Pgas/p, 
and s = (F - 1)"^ ln(Pgas/p^). The GRMHD cooling flux is 
simply Ocooi — Fcom where Fcom is the comoving disc flux, 
as given by Eq. 4. 



D2 Net result of the luminosity calculation 

The flnal goal is to solve for the value of Qcooi that satisfles 
Eq. D2. The physical interpretation of this newly derived 
Qcooi is simply the cooling rate corresponding to a disc with 
the vertical structure given by TLUSTY that is heated up 
according to the GRMHD dissipation proflle. Putting ev- 
erything together, substituting Eq. D5 for Qhoat into Eq. 
D2 gives Qcooi as: 



A Sim 
^adv 



Qa 



(D7) 



QaJ^ and Qco™i are computed from Eq. D6 and Eq. 4 re- 
spectively. Qadv uses the TLUSTY vertical structure, and 
is evaluated via Eq. D3. However, the process of locating the 
correct value for Qcooi that solves Eq. D7 is complicated by 
two factors: 

(i) Qadv on the right hand side is a function of Qcooi- 
Qadv indirectly depends on Qcooi through the annuli 
vertical structure (since Tcff ~ [Qcooi/o'sb]^'^* is an annuli 
parameter). The strategy that we employ to flnd the correct 
Qcooi is a bisection method; we start with the two initial 
guesses for Qcooi that bracket the relation in Eq. D7 (which 
we have empirically found to be monotonic with Qcooi). 
We then bisect on this Qcooi interval until Eq. D7 is satisfied. 



/■7r 
9=0 P\ 



(D4) 



where u^i^ represents the pointwise radial velocities mea- 
sured from GRMHD simulations. 

Our goal is to find the value of Qcooi that satisfies Eg. 
D2, given that we can measure Qh™t from the simulations 
and calculate Qadv from Eq. D3. 



Dl Obtaining Q^™t (the GRMHD dissipation 
profile) 

Unfortunately, the GRMHD simulations that we ran were 
not set up to directly keep track of the numerical dissipa- 
tion Qh™t('')- Despite this lack of information, we can still 
estimate the dissipation indirectly by running the argument 
in Eq. D2 backwards; we solve for the vertically integrated 
dissipative heating rate 



^sim 
/heat 



Tadv + ' 



^slm 
^cool • 



(D5) 



(ii) To compute the ds/dr term, we need to know Qcooi 
(or equivalently Tcs) for two neighboring annuli. However 
the process outlined above (in point 1) only lets us solve 
Qcooi for a single annulus. The resolution to this problem is 
to choose a value of Qcooi for the outermost annulus as a 
boundary condition. If there are a total of A'^ annuli, then 
given Qcooi for the A'^*'' annuli, we can solve Qcooi for the 
{N - 1)*'^ annuli. This process can be iterated to obtain 
Qcooi for all remaining interior annuli. We set Qcooi ~ Qhcat 
as the boundary condition for our outermost annuli since 
far into the disc, energy advection becomes negligible (in 
other words, Qadv — far out in the disc, which implies 
Qcooi -> Qhoat from Eq. D2). 



To maintain consistency with the simulation cooling flux Q^"") 
(as calculated by Eqs. 3 and 4, which only considers bound disc 
fluid) , the integral in Eq. D6 is likewise only taken over the bound 
fluid. 
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Running through the above two steps allows us to 
solve for Qcooi, and hence obtain a second measure of the 
disc luminosity (labelled as GRMHD2 in all plots and ta- 
bles). Note that the GRMHD2 model is more self-consistent 
with TLUSTY in Fig. 14. The TLUSTY2 model (computed 
from the newly derived GRMHD2 disc luminosities) and the 
intrinsic GRMHD entropy profiles agree fairly well when 
r > nsco. 



